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ABSTRACT 

We have investigated the planetesimal accretion rate onto giant planets that 
are growing through gas accretion, using numerical simulations and analytical 
arguments. We derived the condition for gap opening in the planetesimal disk, 
which is determined by a competition between the expansion of the planet's Hill 
radius due to the planet growth and the damping of planetesimal eccentricity due 
to gas drag. We also derived the semi-analytical formula for the planetesimal ac- 
cretion rate as a function of ratios of the rates of the Hill radius expansion, the 
damping, and planetesimal scattering by the planet. The predicted low planetesi- 
mal accretion rate due to gap opening in early gas accretion stages quantitatively 
shows that "phase 2," which is a long slow gas accretion phase before onset of 
runaway gas accretion, is not likely to occur. In late stages, rapid Hill radius 
expansion fills the gap, resulting in significant planetesimal accretion, which is as 
large as several M e for Jupiter and Saturn. The efficient onset of runaway gas 
accretion and the late pollution may reconcile the ubiquity of extrasolar giant 
planets with metal-rich envelopes of Jupiter and Saturn inferred from interior 
structure models. These formulae will give deep insights into formation of extra- 
solar gas giants and the diversity in metallicity of transiting gas giants. 

Subject headings: planetary systems: formation - solar system: formation 



1. Introduction 



Models of the interior structure of Jovian planets in our solar system suggest that 
Jupiter and Saturn would contain much more amount of heavy ele ments in their envelopes 
than that assuming the solar metallicity (ISaumon fc Guillotl I2004J ). This may imply that 
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significant amount of planetesimals was accreted onto the planets togeth er with gas accretion 



from the protopla netary disk. However, the orbital calculations (e.g.. iTanaka fc Idalll997l ; 
Zhou &: Lirul2007l ) showed that the coupled effect of excitation of planetesimals' eccentricities 
due to gravitational scattering by the planet and their damping by aerodynamical and/or 
dynamical drag tends to open up a gap in the planetesimal disk, resulting in truncation of 
planetesimal infall onto the planet's gas envelope. 



Pollack et al.l (119961 ) a priori assumed the maximally efficient planetesimal accretion 



during gas accretion phase. As a result of increase in the planet's mass, the width of its 
feeding zone, which is proportional to cubic root of the mass, expands. They assumed that 
planetesimals in the expanded zones are accreted with the fastest rate for circular orbits of 
planetesimals. Their assumption can be consistent with the anticipated metal-rich envelopes 
of Jupiter and Saturn, but is inconsistent with the eccentricity excitation and gap formation 
shown by the above orbital integrations. 

Furthermore, the assumption of the maximal planetesimal accretion results in long 
"phase 2" that is very inefficient gas accretion phase before onset of runaway gas accretion. 
As explained in §2, envelope contraction starts when core's mass (M c ) becomes larger than a 
critical core mass (M C5 h y dro)- For M c > M C; h y dro, pressure gr adient no more supports envelop e 
gas hydrodyn a mical ly against the increased core's gravity (lMizunolll980l ; llkoma et al.ll2000l ). 
Pollack et al.l (119961 ) showed that heat generation due to the assumed planetesimal accre- 
tion associated with gas accretion supports the envelope quasi-hydrodynamically (in other 
words, it increases the critical core mass; eq. pQ), after the onset of envelope contraction. 
The quasi-hydrodynamical state is called "phase 2" and it may last for more than Myrs. 
However, the in efficient gas ac cretion would be inconsistent with the ubiquity of extrasolar 



giant planets (llda fc Linl 120081 ) . 



Recently, likelihood of phase 2 is re-addressed. iFortier et al.l (120071 ) showed that even 
if the gap form ation is neglected, m o re realistic planetesimal acc retion rate based on oli- 



garchic growth dKokubo fc Ida 



1998 



2002J; iThommes et al.l 120031 ) significantly suppresses 



the duration of phase 2. IZhou &: Linl (120071 ) showed that planetesimals around a protoplanet 



with M c ~ M c hydro are gravitationally shepherded and cannot be accreted. It suggests non- 
existence of phase 2. They also showed that the anticipated metal-rich envelopes of Jupiter 
and Saturn is not inconsistent with it, because such shepherding occurs only in early stages. 
The planet starts accreting planetesimals when its mass becomes comparable to that of gas 
giants because the planetesimals trapped in mean-motion resonances are released by reso- 
nance overlapping due to the planet's mass increase. The planetesimal accretion no more 
halts gas accretion onto such a massive planet. Through numerical simulations with two 
different simple gas accretion prescriptions, they estimated that total accreted mass can be 



- 3- 



as large as several earth masses. 



The idea by IZhou fc Lid (120071 ) reconciled the efficient formation of gas giants with 
the anticipated metal-rich envelope. However, since they showed only numerical results 
with limited prescriptions for gas accretion, it is not clear that the accreted planetesimal 
mass for more realistic gas accretion rate is as much as that they obtained. Furthermore, 
they discussed suppression of phase 2 only qualitatively. As shown below, the planetesimal 
accretion rate does depend on gas accretion speed as well as the planet's mass. Hence, 
for evaluation of amount of planetesimal infall for more realistic gas accretion models and 
quantitative discussion on the suppression of "phase 2," general formulae for gas accretion 
rate as a function of the planet's mass (M) and its increase rate (M) is needed. 

In the present paper, through orbital integrations, we clarify the physical mechanism 
to determine the planetesimal accretion rate and derive detailed semi-analytical formulae 
for the accretion rate as a function of M and M. We find that the total infall mass of 
planetesimals can be as much as several earth masses even for more realistic gas accretion 
models and quantitatively show that phase 2 is unlikely to occur. 

The outline of this paper is as follows. We summarize gas accretion processes onto plan- 
ets in §2. The method of our calculation and initial setup is described in §3. With artificial 
simple gas accretion models, we clarify intrinsic physics that determines the planetesimals 
accretion rate and derive semi-analytical formulae for the accretion rate (§4.1 to 4.3). Ap- 
plying the formulae to realistic gas accretion models, we discuss the metallicity of Jupiter 
and Saturn envelope (§4.4). We also discuss "phase 2" and find that phase 2 is not likely to 
occur (§4.5). The conclusion is in §5. 



2. Gas Accretion Onto a Core 

As mentioned in §1, when the core mass becomes larger than a critical core mass, pres- 
sure gradient no more supports envelope gas hydrodynamically against the core's gravity 
and hydrostatic atmosphere does not exist. After that, heat generation due to gas envelope 
contraction itself supports the envelope against dynamical collapse, thus the envelope un- 
dergoes quasi-static contraction. The contraction allows gas inflow from the disk into Bondi 
radius of the planet, so that the contraction rate is almost equal to gas accretion rate of the 
planet. 

Here, we briefly summarize the prescriptions for this process for later use. The critical 
core mass depends on planetesimal accretion rate onto the core (M c ) a nd the grain opacity 



K gI ) associated with the disk gas. Based on a series of numerical models, llkoma et al.l (120001 ) 
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found that the critical core mass for break-down of hydrostatic atmosphere is 

0.2-0.3 



M, 



c, hydro 



10 



lO^Meyr- 1 



0.2-0.3 



where 



lcm 2 g _1 ) is the grain opacity given by lPollack et al.l (119851 ). who assumed dust 
grains with interstellar abundance and size distributions. Faster accretion and higher opacity 
(relatively large M c and K gr ) result i n a warmer planetary atmosphere and an enhanced 
pressure gradient, so M c hydro is larger (lStevensonlll982l ; llkoma et al.ll2000i ). 



Pollack et al.l (119961 ) assumed the most efficient planetesimal accretion induced by ex- 
pansion of feeding zone due to increase of the planet mass, the rate of which is ~ 10 _6 A / / e yr _1 , 
for M c ~ 10M e . When the planet with M c ~ 10M e becomes isolated consuming planetesi- 
mals in its feeding zone, gas envelope starts contraction and the induced planetesimal accre- 
tion from the expanded region of the feeding zone increases M C) h y dro up to ~ M c (eq. p]) and 
stall gas accretion. This self-regulated process works on more than Mrys until M c exceeds 
~ 20M e . This is called "phase 2." However, as we show in §4.5, the rate of the planetesimal 
accretion induced by gas accretion is not generally large enough to maintain phase 2. Then, 
gas accretion dominant phase starts. 

For M c ~ M Ci h y dro, heat generation due to planetesimal accretion marginally equili- 
brates with the core's gravity. In the quasi-static contraction stage, heat generation due to 
envelope contraction marginally equilibrates with the gravity of the planet with total mass 
M (including envelope mass). The Kelvin- Helmholtz contraction timescale is equivalent to 
planet mass increase timescale r g)acc = M/M. Replacing M c and M c by M and M/r gjacc in 
eq. flD, 7V acc is given by 



10 



7 / 



M 



10M 



-(2.3-4) 



Ken- 



(2) 



Deta i led numerical simulatio ns of quasi-static evolution of the gaseous envelope (llkoma et al. 
2000l ; llkoma &: Genda 1 120061 ) show consistent results at the onset of runaway ga s accretion 
in which the envelope and core masses are nearly equal. Although IPodolakl (120031 ) suggested 
K gT ~ 0.0lKg r through the numerical simulations of coagulation and sedimentation of dust 
grains in the atmosphere, the amount and size distr ibution of dust gr a ins in the atmosphere 



are highly uncertain. Here, we adopt the results by llkoma fc Genda I (120061 ) with n gr — ,,„, . 

M 



T, 



g,acc 



10 
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V ioM ff 



-3.5 



yrs. 



(3) 



fiducial "realistic" gas accretion model. 
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When eq. (j3j) is extrapolated to large M (> 100 M©), it may give un realistically fast sup- 



ply of gas from the disk. Hence, we limit the gas accretion rate as bellow. iTanigawa &: Watanabe 



( 120021 ) showed through two-dimensional local hydrodynamic simulations, the mass infall to 



the circumplanetary subdisk from the protoplanetary disk is limited by 



M *. f a / M 

— ~ 6 x 10" 4 
M 



, / a \-i-s / M \ 1 



where / s is a scaling factor for disk gas surface density defined by eq. 0. We use this limit 
with f g = 0.7. 

Another limit is Bondi gas accretion, the rate of which is given by M = irr B p gas c s , 
where p gas = T, g /(2H) is the spatial density of gas disk and H is the disk scale hight, 
tb = 2GM/c 2 s is the Bondi radius, c s = HQ is the sound speed and Q = a/ GM & / a 3 is 
the Keplerian angular velocity. Adopting the temperature distribution in the limit of an 
optically thin disk (Hayashi 1981), T = 2.8 x 10 2 (r/lAU) _1/2 K, the Bondi gas accretion 
limit is 

0-7 X10- 3 ( — ) — yr" 1 . (5) 



M V5AU/ \M ( _ 

As figure [1] shows, the timescale (= M/M) in eq. (j4j) is generally longer than the Bondi 
accretion timescale, so an actual lower limit for gas accretion timescale is given by eq. (jl]). 



3. Calculation Setup 

3.1. Orbital Integration 

We numerically calculate the orbital evolution of a swarm of planetesimals in the vicinity 
of a protoplanet's orbit embedded in a gaseous disk. The protoplanet grows accreting gas 
with a given rate. The planetesimals are treated as massless test particles and neglect their 
interactions. We assume that the protoplanet has a fixed circular orbit. 

The planetesimals' orbits are affected by the gravitational force of the growing proto- 
planet and drag force from disk gas, 

f = _ V ~ V S™ (Q) 
J gas 5 \ u ) 

1 damp 



where v and f gas are the velocity of a planetesimal and disk gas. The gas motion is a 
circular Keplerian motion. In some runs, we adopted sl ightly slower rotati on speed of the 
disk gas due to radial pressure gradient in disk gas (e.g., lAdachi et al.lll976l ). which induces 
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inward migration of planetesimal orbits. However, we found that the inward migration hardly 
changed the results of planetesimal accretion rate onto the protoplanet. We here show the 
results without the inward migration. We set a damping timescale of the gas drag, r damp 
(= e/e), as a constant parameter for all the planetesimals throughout a run in order to make 
clear the effect of the damping force. 



We follow the prescription of gas surface density distribution by llda fc Linl (12004 ). 



2l0f 9 



5AU 



-3/2 



gem 2 , 



(7) 



where f g is the scaling parameter and f g = 0.7 corresponds to the gas surface density of the 
minimum mass solar nebular model, £ 9i mmsn (jHayashilll98l[ ). For simplicity, we neglect a 
gap i n the gas disk, which may be opened up by the perturbations from a massive protoplanet 
(e.g.. lLin &: Papaloizoulll993l ) and assume the above unperturbed S 9 eve rywhere. With this 



S ff , a given value of Td amp corresponds to individual planetesimal mass (lAdachi et al.l 11976 
Tanaka fc IdaHl999h . 



m = 3 x 10 17 f 3 — 
Ja V0.1 
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a 
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where a, e, and p p \ are semi-major axis, e ccentricity, and materi al density of planetesimals, 
respectively. If gravitational drag (e.g., iTanaka fc Wardl 12004 ) is considered in stead of 
aerodynamical gas drag, 



m 



4.5 x 10 26 /, 



-1 / 7~damp \ 



a \ ■ 
5Mj) 



(9) 



although in this case, interactions among the planetesimals could be important. 



The orbi ts of planetesimals are n umerically integrated by using t he fourth-orde r Her- 
mite scheme (jMakino Sz Aarseth 1992 ) with the hierarchical timestep (Making 1991 ). The 
equation of motion of particle k is given by 



d 2 r k 
dt 2 



-GMr. 



r k GM(r k - r) GMr 



0] 



\r k 



+ f 



gas 5 



(10) 



where M and r is the mass and position of the protoplanet. The first term to the last 
one represent the gravity from the central star, the gravitational perturbation from the 
protoplanet, the indirect term and the gas drag force, respectively. We set M* = M . 



When a planetesimal contacts the surface of the protoplanet, the planetesimal is removed 
after recording the collision. The planet mass is unchanged. The physical radius of a 
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protoplanet is determined by its mass and internal density p as 



R 



( 




(11) 



We set p = lgcm 3 in all simulations. The dependence of the planetesimal accretion rates 



on p will be discussed in §4.4. 

Although we neglect gravitational forces of planetesimals, mass of planetesimals is spec- 
ified in order to calculate the amount of mass accretion onto the protoplanet (regarding 
"effective" mass for gas drag force, see below). Assuming equal- mass planetesimals, they are 
initially distributed in the range a in < a < a out to satisfy the surface mass density 



where fd is a scaling factor. As is the case for f g , fd = 0.7 corresponds to MMSN. The 
inner and outer boundaries are a in = a p (l — 5/tf) and a out = a p (l + 10/if), where a p is the 
semi-major axis of the protoplanet and h{ is the reduced Hill radius for final mass of the 
planet (M f ). The reduced Hill radius of a protoplanet h is Hill radius r H scaled by a p , 



In all numerical simulations, we adopt a = 5AU, Mf = Mj (Jupiter mass), and fd = 
f g . Accordingly, h { = 6.8 x 1CT 2 , a in = 3.3 AU, and a out = 8.4AU. We will derive the 
dependences on a, Mf, and fd(= f g ) by analytical arguments and discuss the results with 
different parameter values. Total mass of planetesimals within the region a m < a < a ou t is 
~ 20 /dM e . The number of planetesimals in most runs is iV = 20000. With an assumption 
that planetesimals have an equal mass, their individual mass is m p i ~ 1.1 x lO^/^M^ = 
6.6 x 10 24 /d g. In our simulations, we specify r damp = 10 6 , 10 5 , 10 4 yrs and oo (gas-free case), 
independent of the values of m p i. Except for the gas-free case, the above values of m p i is much 
larger than the values in eq. (jSJ) for the given Td am p, so planetesimals that we use correspond 
to "super particles" representing many smaller planetesimals. Since we neglect interactions 
among planetesimals, such "super particles" treatment is not inconsistent. Initial eccentricity 
and inclination of planetesimals are taken as eo = io = 0.001 for all simulations. Since e and 
i are quickly pumped up by perturbations from the protoplanet, the choice of eo and io does 
not affect results. 




(12) 




(13) 
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3.2. Growth of a Protoplanet 

Since we consider the phase after isolation of protoplanets, we assume that the growth of 
the protoplanets is dominated by accretion of surrounding disk gas but not by planetesimals. 
As we will show later, this assumption is valid, because amount of accreted gas is much larger 
than the anticipated amount of accreted planetesimals. 

In §2, we described the prescription for gas accretion. In the numerical simulations, 
we use simple artificial gas accretion models in order to make clear what conditions regu- 
late the planetesimal accretion rate. From the results with the artificial models, we derive 
semi-analytical formulae for the planetesimal accretion in general forms (§4.3). Applying 
the formulae to the more realistic gas accretion rate in §2, we calculate the total mass of 
planetesimal infall into the envelope of Jupiter and Saturn in §4.4. 

The simple artificial gas accretion models are expressed by 

dM =aM^ (14) 



dt 



where a is the integration constant determined by the bound ary condition . We set the 



condition as M = M for t = and M = M f for t = t f . Following IZhou fc Lid (120071 ). we set 
the protoplanet at 5AU with its initial mass Mo = 5.67M© and final mass Mf = Mj (Jupiter 
mass). We adopt tf = 10 5 yr for numerical simulation, following their nominal cases. Th e 



growth with p = 2 and correspond to the Bondi and linear models in IZhou fc Lid ( 120071 ). 
Figure [2] shows the evolution of the mass of a protoplanet by accreting gas for p = 2, 1, 0, —2. 
Here we assume that Mq > M C) h y dro- The consistency of the assumption is checked in §4.5. 



4. Results of Orbital Calculation 

4.1. Overall Evolution 

Figure [3] shows the snapshots of the distributions of planetesimals on the b-e/h plane, 
where h is defined by eq. (fT31) . The scaled orbital separation b is defined by 

b= a —^, (15) 
na p 

where a p is semimajor axis of the protoplanet. The protoplanet is fixed at the origin (i.e., 
b = and e = 0). The growth rate of the protoplanet M is oc M 2 . To avoid busy plots, 
we show only 1000 planetesimals in this figure. We integrate the evolution of planetesimals 
for 3 x 10 5 yrs. Since tf = 10 5 yrs, we set M = for t = 1-3 x 10 5 yrs, which corresponds to 
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termination of gas accretion due to gap formation in the gas disk, although we neglect the 
effect of gas density depletion on drag force. The damping timescale is Td am p = 10 4 yr. 

In figure [31 we also drew the Jacobi energy Ej, 

Ej = l -{{e/hf + (z/h) 2 ) - h 2 + 9 - + 0(h). (16) 

In the figure, we also include higher order terms of h in Ej. In the circular restricted three- 
body problem, Ej is conserved between before and after scattering by the protoplanet (on 
average, both (e 2 + i 2 ) and b 2 incre ase) . Since only plan etesimals with Ej > can enter the 



Hill sphere of the protoplanet (e.g.. lHayashi et al.lll977l ). we regard the region Ej > as the 



feeding zone of the protoplanet. When e/h, i/h < 1, the width of feeding zone is 

In the top panel (t = lOOOyrs), the planetesimals in the vicinity of the protoplanet are 
scattered and their e and b increase along a constant Ej curve. The planetesimal eccentric- 
ities e are damped in the panel of t = 3 x 10 4 yrs because of r damp = 10 4 yrs. Since the gas 
drag damps e keeping b almost constant, all the planetesimals except for those trapped in 
horseshoe orbits go out of the feeding zone. As the protoplanet grows up, its feeding zone 
expands. Since b oc M -1 / 3 , b of planetesimals decreases with time, but scattering opposes 
it. Since M oc M 2 , the expansion accelerates with time. Eventually, the expansion over- 
whelms the gap opening due to a coupling effect of scattering and gas drag damping, so that 
planetesimals go into the feeding zone in the panel of t — 1 x 10 5 yrs. After t — 1 x 10 5 yrs, 
the increase of M stopped, so a gap in the planetesimal disk is again produced (the bottom 
panel). Planetesimals that have sufficiently large b are captured in proper mean motion res- 
onances. Although we neglect inward migration due to slightly slower rotation of gas than 
Keplerian rotation due to pressure gradient, damping of e results in small decrease in b due 
to angular momentum conservation. Such inward migratio n causes resonance trapping. The 



evolution of the gap width is consistent with the result by lZhou fc Lin! (120071 ) . 



Thus, after initial relaxation, planetesimals are shepherded and planetesi mal accretion 



rate i s very low, until efficient planetesimal accretion re-starts in late stage. IZhou fc Lin 



( 120071 ) showed through orbital simulation of planetesimals that planetesimal accretion oc- 
curs only in late stage of gas accretion. They suggested that in late stage, planet's mass 
becomes large and the mean motion resonances overlap to release planetesimals captured in 
the resonances. So, they concluded that mass of the protoplanet controls the accretion rate 
of planetesimals onto the protoplanet. This effect indeed determines supply of planetesimals 
to the regions near the feeding zone. However, whether the planetesimals near the feeding 
zone are shepherded or not is determined by gap opening that is a result of competing pro- 
cesses of the feeding zone expansion and scattering/damping, so the values of M play an 
important role as well as M. 



-10- 



4.2. Dependence of Planetesimal Accretion Rate on Planet's Mass 

The evolution of the planetesimal accretion rate for rd am p = 10 6 , 10 5 , 10 4 yrs and gas- 
free case is plotted as a function of M in Figure HI The four lines in each panel represent 
various gas accretion models (p = 2, 1,0, —2). Initial mass of the protoplanet M is set as 
5.67 M ffi . Starting with 20,000 planetesimals (i.e., planetesimal mass ~ 1.1 x 10~ 3 /dM e ), we 
calculated for 10 5 yrs (the growth timescale if = 10 5 yrs). 

To see the dependence on M more clearly, we plot the scaled planetesimal accretion 
rate M/R 2 , where R is the physical radius of the protoplanet. Through the numerical 
simulations, we found that planetesimals are likely to experience 2-D accretion rather than 
3-D. The 2-dimensional accretion rate is 

dM orv f v «°\ / 32vrGp ^ 2 

~df ~ 2iffid t J " rcl = v ~^- EdR ' (17) 

where S^, t> esc = ^/2GM/R and v Ie \ are the surface density of the planetesimals, escape 
velocity from the protoplanet 's surface and relative velocity between the protoplanet and 
planetesimals, respectively. The scaled accretion rate (M/R 2 ) is determined by effective T, d 
in the feeding zone for fixed internal density of the protoplanet (p). In our simulations, the 
total mass of the planetesimals is not significantly decreased, so the effective is determined 
by scattering by the planet, gas drag, and Hill radius expansion due to the planet growth. 

Figure H] shows that the planetesimal accretion rates for fd = 0.7. The planetesimal 
accretion rates as a function of M depends on the parameter p. For p = —2 and 0, the scaled 
planetesimal accretion rate decreases with M, which suggests that a gap in the planetesimal 
disk is formed when M becomes large. On the other hand, for p = 2, the protoplanet may 
grow so fast in the late stage that the feeding zone expansion overwhelms the gap formation, 
as shown in figure [3j The dependence on p implies that the accretion rate is not a function 
solely of M, but depends on M as well as M, because p determines the M-M relation. 



4.3. Dependence of Planetesimal Accretion Rate on Gap Formation 

Parameters 

Here we show that competition among the feeding zone expansion, scattering, and 
eccentricity damping regulates flux of planetesimals across the boundaries of the feeding 
zone (Ej = 0), that is, the planetesimal accretion rate. We consider change rates of b 2 and 
(e/h) 2 of planetesimals (we neglect the contribution from i because i is usually correlated to 
e and i < e). 
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Evolution of planetesimals on the b 2 -(e/h) 2 space due to gravitational scattering by 
the protoplanet, damping of eccentricity by gas drag, and expansion of Hill radius by mass 
increase of the protoplanet is expressed by the change rates, f sca t, ^dam P , and t>n, on the 
plane. Since 6 oc h^ 1 oc M -1 / 3 , 



dtf 

~dt 



(growth) = (bh) 



,dh 



dt 



2 l9 M 
-b z — 
3 M 



T, 



g,acc 



where r„ 



= M/M is the timescale of planet mass increase. In the last equation, we 
used b ~ 2\^3, which is the location of the feeding zone for e/h < 1, for simplicity. With 

''"damp ^/^j 



ld(e/hy 



amp 



(damping) 



(e/hf 



(19) 



2 dt 7~damp 

The factor (1/2) in the definition is added for the more simple form of the final expression 
and better fit with numerical results. Evolution due to the scattering is increase of b 2 and 
(e/h) 2 on average, along a constant Ej curve (Ej ~ 0). The corresponding change rate is 



db 2 , 

Wscat = —(scattering) 
dt 



Ad(e/hf 
3 df~ 



(scattering). 



(20) 



Assuming long-range gravitational interaction with (e/h) < 1, linear calculation (IGoldreich fc Tremaine 
1982 : Hasegawa fc NakazawalE99ol ) showed that b of a planetesimal is increases by 5b ~ 306~ 5 
during each encounter. Numerical calculation showed that for b ~ 3-4, 5b is overesti- 
mated by a factor ~ 10 (lldal Il990l ). Since the scattering occurs at every synodic time 



T ~ 

- 1 - syn — 



27ia p /(lbr n Q K ), 



,db 

^scat = 26— (scattering) 
dt 



2b 



0.15b 



T, 



syn 



6 (3/2)6/i h 
— y 1 ' ~ 0.22 — . 
b 4 T K T k 



(21) 



where Tk = 27r/f2K is Keplerian period and 6 ~ 2^3 is again used. 



Since the feeding zone is determined by the values of Ej and the scattering does not 
change the values, the condition of gap opening would be fdamp ^ ^h- If inward migrations of 
planetesimals due to gas drag or type I migration of the protoplanet is considered but growth 
of a protoplanet is neglecte d, the gap formation condition is similarly derived, replacing 
by db 2 /dt due to gas drag (ITanaka &: Idalll997l ) or type I migration (ITanaka fc Idalll999l ). 
These effects can suppress g rowth of the protoplanets before they attain their isolation masses 
jTanaka fc Idalll997l . \l99% . 



When fdamp ^ ^h, the gap is not created and planetesimals are engulfed by the expand- 
ing feeding zone. The engulfment rate would be determined by vu/v SCSLt , because vu and 
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/ v sca t have opposite directions to each other in the b 2 components. Thus, it is expected that 
for Vdamp < wh, the accretion rate would be regulated by 

e = i J* i ~ srh-^ ~ 4 1 (^) 3/2 f^LY 1/3 f^Y 1 (22) 

while for t>dam P > ^h> the accretion rate would be regulated by 

„ = I V * I ~ 8 T damp ^ q g f ^damp \ ^ / T g , aC c \ / M\ ^ /^\ 3 / 4 f , 

' Vamp' {e/hf r gi acc ' V10*yra/ VIOVs/ VM®/ ^5AU7 ' V ; 

where we used eq. (13131) in Appendix for (e/h) 2 . Since r gacc = M/M and Td am p and /i are 
functions of M, the planetesimal accretion rate would depend on M as well as M. We 
show that the numerical results agree with the above argument and derive formulae for the 
planetesimal accretion rate as a function of £ and rj. 

Figure shows the evolution of the scaled planetesimal accretion rate as a function of 
£ = vu/v scllt for p = 2, 1, and —2 in the cases of Td am p = 10 6 , 10 5 , 10 4 yrs and gas-free case. 
As suggested in the above discussion, figure [5] shows that in the ranges of i] > 1 [equivalent ly, 
£ > 5(r dam p/10 5 yr) _1/2 (M/M e )~ 1/6 (a p /5AU) _3/4 ], the scaled planetesimal accretion rate is 
independent of planet gas accretion models with different p (different M—M relations) and 
different Td amp . This confirms that the accretion rate is determined by £ for rj > 1 (non-gap 
cases). The planetesimal accretion rate in this case is given by 

with a ~ 0.8 and j3 ~ —6 that are obtained from our numerical results by the least square 
fitting. The fitting line, eq. (1241) . is expressed by thick solid lines in the plots. For r\ < 1, the 
accretion rate declines, which corresponds to gap opening in the planetesimal disk. 

Figure [6] shows the evolution of the scaled planetesimal accretion rate as a function of 
V = fH/^damp- In the range of r] < 1, the scaled planetesimal accretion rate is independent of 
planet accretion models with different M-M relations and different Tdamp- This confirms that 
the accretion rate is determined by 77 for 77 < 1. From our numerical results, the planetesimal 
accretion rate in this case is given by 



dM sohd 

dt 

with a ~ 1.4 and j3 — —6. 



W#V/-(— y^ey*- 1 . (25) 

\ -ft® / V ^damp / 
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When the planet's mass has grown to M, the total mass of planetesimals that infall in 
the envelope (M soiid (M)) is obtained by integrating 



dM. 



solid 



dM 



^-^solid ^g,acc 

~dt W 



(26) 



from to M. In figure U\ M SO Hd(M) evaluated by the above semi-analytical formulae is 
compared with that obtained by orbital calculations for individual gas accretion models in 
cases of Td am p = 10 6 yrs, 10 5 yrs and 10 4 yrs. The semi-analytical formulae well reproduce the 
results of orbital calculations except for early stages in which M so iid is so small that statistical 
fluctuation i s larg e. The formulae also reproduce numerical results in their figure 9a in 
Zhou fc Linl J2007h . 



4.4. Application to Jupiter and Saturn 

In the preceding subsection, we investigated planetesimal accretion onto growing pro- 
toplanets with artificial gas accretion models and obtained semi-empirical formulae of the 
planetesimal accretion rate. Applying this formulae to the more realistic gas accretion models 
in §2, we discuss the metallicity of envelopes of Jupiter and Saturn. 

Integrating eq. fl26l) with eq. ([3]) to Mf, we estimate total mass of the accreted plan- 
etesimals in the cases of Jupiter (Mf = 318M®, a, p = 5.2 AU) and Saturn (Mf = 95M e , 
a = 9.55 AU). The evolution of M solid is plotted in figure [SJ The three curves show the 
results with Td am p = 10 4 , 10 5 and 10 6 yrs. It is likely that gas giants were inflated during 
gas accretion phase. For a fixed M, dM/dt oc fd^fpR 2 oc fdVR (eq. [IZ])- In the figure, 
we plot the accreted planetesimal mass M* olid for R = 2R 1 and fd = 2, where Ri is the 
physical radius for mass M and p = lgcm -3 . For other R and fd, the accreted mass is 
M solid = (i?/2i? 1 ) 1 / 2 (/,/2)M s ; iid . 

All the results show similar qualitative features of evolution of M so i i( j. Planetesimal 
accretion is inhibited in early stages by gap formation, but rapid planetary growth due to 
gas accretion in later stages allows planetesimal accretion. With Td am p = 10 6 yrs, M so nd — 
6(-R/2i?x) 1//2 (/ d /2)M e both for Jupiter and Saturn. For shorter r damp , M so i id is smaller due 
to easier gap formation. We also did calculations starting from different core masses. The 
resultant M so iid hardly changed, because dM so na/ dM is negligibly small when M is small 
and gap is opened. The am ount of predicted M sr ,iid ca n be as large as that inferred from 



the internal structure model ISaumon fc Guillotl (120041 ). if the planets are inflated and/or 
relatively large fd is considered. 



For the same M, p and fd-, M so iid is larger for larger a p . Although Saturnian mass is 
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1/3 of Jovian mass, our model predicts that the mass of planetesimals falling into Saturnian 
envelope is comparable to that into Jovian envelope. More detailed internal structure models 
will test our prediction. 



4.5. Phase 2 



So far, we have assumed that gas accretion immediately starts when M c exceeds M c hydro 
without undergoing "phase 2." In the previous subsection, we predicted the planetesimal 
accretion rate as a function of planetary mass based on the realistic gas accretion model. 
With this accretion rate, we show that "phase 2" is not likely to occur. 



In the nominal model (Jl model) in iPollack et al.l (119961 ). fd — 2.5, a p = 5.2AU and 
M c ~ IOMq. Then, they found that M c ~ 10 _6 M e /yr is maintained during "phase 2" with 
their maximally efficient planetesimal accretion model. As shown in eq. ([IT) , this M c can 
marginally support gas envelope around a 10M e core. 

First, we derive the condition for gap opening with a realistic r gjacc given by eq. ([3]). 
Substituting eq. (j3J) into eq. (1231 . 



7] ~ 0.8 x 10 



—6 / ^damp \ 

V 10 4 yrs J 



1/2 



M 

Mm 



3.3 



5AU 



3/4 



(27) 



With a- = 5.2AU, M ~ M c ~ 10M e and r d 



amp 



10 6 yrs, we obtain ?; ~ 2 x 10 2 C 1. 



Then, the gap should be opened up. Our formula for rj < 1 gives 



M solid ~2.2xl0- 6 ^ 



Pp 



lgcm 



-1/6 



\ 7/10 
Tdamp \ 



( T g,acc \ 
10 4 yrs / 



-7/5 



M 

Mm 



13/30 



( a P 

V5AU 



21/20 



M e /yr. 
(28) 



Substituting eq. into this equation, 

-1/6 



solid 



0.9 x 10" 14 /d 



lgcm 



-3 



\ 7/10 
Tflamp \ 



10 4 yrs / 



M 



16/3 



f a P \ 

V5AU/ 



21/20 



Me/yr. (29) 



For r damp = 10 6 yrs, f d ~ 2.5, a p = 5.2AU and M ~ 10Af©, M r ~ 1.1 x lQ- 7 M m /yr, which 
is one order smaller than the planetesimal accretion rate that IPollack et al.l (119961 ) assumed. 



We examine the possibility of phase 2 for other fd and a p . For phase 2 to occur, M c 
must be maintained to be as large as M for M c ~ M r > , ydl . . Core mass can be approximately 
ident ified by core isolation mass beyond the ice line (IKokubo fc Idal 1 19981 . |2002| ; llda fc Lin 
2004T ). 



lvl CASO 



3/4 



(30) 
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From eq. ([I]) with the exponent derived by assuming eq. ([3]), the accretion rate required by 
occurrence of phase 2 is 



M 



solid,2 



— j AWy r ~3xio-/r 5 (^)"M«/y r . (3D 

Substituting M Cjiso into M in eq. (j2"Tj) . 

^""^fer^r- (32> 

So, 77 < 1 is equivalent to 

For this range of fd and a p > 3AU (the ice line), eq. fl28l) with M replaced by M c iso is always 
smaller than M so iid,2 given by eq. fl3Tl) (see figure [9]). For 77 > 1, on the other hand, 

M. olld = 1.5 x 10-A f^V" 6 f^V 4/5 f^V 2 ' 5 m" /5 M e /y, (34) 
Jd \lgcm- 3 J Vl0 4 yrs/ V M e/ V5AU/ ® / ' y V ; 

In the range of fd and a p that satisfy 77 > 1, eq. flM|) can reach M S oii H.2 only at a p > 15AU 
and fd ~ 1, in which gas giant formation is unlikely (llda fc Lin! 12004 ) . Thus, the predicted 
M c never reaches the values required for phase 2. We conclude that phase 2 is not likely 
to occur for formation of giant planets. This conclusion is consistent with the ubiquity of 
extrasolar gas giant planets. 



5. Conclusion 

We have investigated the planetesimal accretion rate onto growing giant planets through 
numerical simulations and analytical arguments. The planet mass (M) is increased with 
assumed gas accretion rate onto the planet, and orbits of planetesimals in the vicinity of 
the planet's orbit are integrated with the effect of gas drag, but without self-gravity of the 
planetesimals. 

We first performed simulations with several different artificial gas accretion rates to clar- 
ify intrinsic physics determining the planetesimal accretion rate. A gap in the planetesimal 
disk is opened by a coupling effect of gravitational scattering by the planet and gas drag 
damping. Here, the gap formation means that most planetesimals are get out of the feeding 
zone of the planet. The scattering increases both e and b keeping Jacobi energy constant, 



where e is orbital eccentricity and b is difference in semimajor axis between the planet and 
the planetesimals. Changes in e and bh are of the same order, where h is reduced Hill radius 
defined by (M/3M*) 1 / 3 . Since the gas drag predominantly damps e after the scattering, 
the gap is formed. On the other hand, the width of the feeding zone is proportional to h. 
Thus, the planet growth inhibits gap formation and competes with the scattering/damping 
process. 

We derived the condition for the gap formation by comparison between the eccentricity 
damping rate (v damp) an d the rate of expansion of the feeding zone due to the planet growth 
(^h). When vn/vdam P > 1, the gap is not formed. Then, the planetesimal accretion rate 
(dM so nd/dt) is scaled by the ratio of the scattering rate t> scat to t>H- The numerical results 
are fitted as 



where R is physical radius of the planet and fd is a scaling factor for surface density of 
the planetesimals (eq. [12] ). When the gap is formed (v^/v^mp < 1), the accretion rate is 
significantly depleted. We found that the accretion rate is scaled by fH/^damp as 



Applying these formulae to the more realistic gas accretion models described in §2, we 
found the followings: 

1. In early stages when M ~ O(10)M e , a gap is opened in the planetesimal disk. The 
planetesimal accretion rate is smaller than that required for phase 2 to occur. This 
ensures efficient formation of gas giants, which may be consistent with the ubiquity of 
extrasolar giant planets. 

2. In later stages (M > O(100)M e ), the expansion of the feeding zone overwhelms the 
gap opening process, so the gap is filled. Then, the planetesimal accretion becomes 



3. The amount of infalling planetesimals into the envelopes of Jupiter and Saturn in the 
late stages can be as large as several M®, which may be consistent with interior models 
for these planets. 

In this "realistic" model, we assumed that planetesimals are infinitely supplied. However, 
if the accreted mass is significant, planetesimals distributed in the regions inside isolated 
strong mean motion resonances can be consumed. In that case, release of planetesimals 




(35) 




(36) 



efficient. 
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from the resonance capture by resonance over lapping due to planet mass increase may also 
become a important factor (jZhou fc Lin! 120071 ). 



Guillot et al.1 (120061 ) pointed out the correlation that the amount of solid components 



of extrasolar transiting gas giants increases with metallicity of their host stars that is pro- 
portional to fd- This trend is consistent with our formulae, because dM so \\&/dt oc /<j. As 
this example shows, the analysis here will give deep insights into formation of extrasolar gas 
giants and their diversity. 



This work is supported by JSPS. 



Appendix 

The magnitude of {e/h 2 ) in §4.3 is determined by a balance between damping due 
to the gas drag and excitation due to the planet's perturbations. Since in the non-gap 
case, planetesimals are engulfed by the feeding zone mainly through the parameter range of 
(e/h) < 1, we use eq. ( ]21~j) for definition of the parameter £. However, gap opening is caused 
by damping of relatively high orbital eccentricity, so we use the formula of excitation of 
planetesimal eccentricity due to the protoplanet's perturbations for (e/h) > 1, in evaluating 
(e/h) 2 . Then the scatte ring timescale is given approximately by Chandrasekahr's two-body 
scattering formula (e.g., Stewart &: Ida 2000 ; Ohtsuki et al.ll2002l ). 



' e,scat 



n p 7r(GM/(ev K ) 2 ) 2 ev K In A' 



(37) 



where In A ~ 3 and n p is spatial density of the protoplanet, which is given by inverse of 
volume of the planetesimal disk in the feeding zone, l/(2na p x A\f?>ha p )(eV]<i/VL K ). Then, 

8 v / 3vr(e//i) 4 , ,T K „ vi ( M \~ 1/3 / a p \3/2 

From r e)Scat = r damp , we obtain 



REFERENCES 

Adachi, I., Hayashi, C, & Nakazawa, K. 1976, Prog. Theor. Phys., 56, 1756 



-18- 



Fortier, A., Benvenuto, O. G. & Brunini, A. 2007, A&A, 473, 311. 

Guillot, T., Santos, N. C, Pont, F., Iro, N., Melo, C. & Ribas, I. 2006, A&A, 473, L21. 

Goldreich, P., & Tremaine, S. 1982, Ann. Rev. Astron. Astrophys. 20, 249 

Hasegawa, M. & Nakazawa, K. 1990, Astro. Astrophys, 227, 619 

Hayashi, C, Nakazawa, K., & Adachi, I. 1977, Publ. Astron. Soc. Jpn., 2, 163 

Hayashi, C. 1981, Prog. Theor. Phys. Suppl, 70, 35 

Ida, S. 1990, Icarus, 88, 129 

Ida, S., & Lin, D. N. C. 2004, ApJ, 604, 388 

Ida, S., & Lin, D. N. C. 2008, ApJ, 673, 487 

Ikoma, M., & Genda, H. 2006, ApJ, 650, 1150 

Ikoma, M., Nakazawa, K., & Emori, H. 2000, ApJ, 537, 1013 

Kokubo, E., & Ida, S. 1998, Icarus, 131, 171 

Kokubo, E., & Ida, S. 2002, ApJ, 581, 666 

Lin, D. N. C, & Papaloizou, J. C. B. 1993, in Protostars and Planets III, ed. E. H. Levy 
and J. I. Lunine (Tucson: Univ. of Arizona Press), 749 

Makino, J. 1991, Publ. Astron. Soc. Jpn., 4, 859 

Makino, J., & Aarsetth, S.J. 1992, Publ. Astron. Soc. Jpn., 44, 141 

Mizuno, H. 1980, Prog. Theor. Phys. Suppl., 64, 54 

Nakazawa, K., Ida, S. 1988, Prog. Theor. Phys. Suppl., 96, 167 

Ohtuski, K., Stewart G. R. & Ida, S. 2002, Icarus, 155, 436 

Podolak, M. 2003, Icarus, 165, 428 

Pollack, J. B., McKay, C. P., & Christofferson, B. M. 1985, Icarus, 64, 471 

Pollack, J. B., Hubickyj, O., Bodenheimer, P., Lissauer, J.J., Podolak, M., & Greenzweig, 
Y. 1996, Icarus, 124, 62 

Saumon, D., Guillot, T. 2004, ApJ, 609, 1170 



-19- 

Stevenson, D. J. 1982, P&SS, 30, 755 
Stewart, G. R. & Ida, S. 2000, Icarus, 143, 28 
Tanaka, H., & Ida, S. 1997, Icarus, 125, 302 
Tanaka, H., & Ida, S. 1999, Icarus, 139, 350 
Tanaka, H. & Ward, W. R., 2004, ApJ, 602, 388 
Tanigawa, T., & Watanabe, S. 2002, ApJ, 580, 506 

Thommes, E. W., Duncan, M. J. & Levison, H. F. 2003, Icarus, 161, 431. 
Zhou, J.-L, Lin, D. N. C. 2007, ApJ, 666, 447 



This preprint was prepared with the AAS IATgX macros v5.2. 



-20- 




10° 10 1 10 2 10 3 



M [M 8 ] 



Fig. 1. — Gas accretion timescal es for a planet w i th m ass M at a = 5AU. Solid line is 
an extrapolation o f the model by llkoma fc Genda I (120061 ). Dashed and dotted-dashed lines 
represent limits by Tanigawa &: Watanabel ( 20021 ) and Bondi accretion. 
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Fig. 2. — Evolution of the protoplanet mass according to the simple power-law gas accretion 
models (M oc M p ). Initial and final masses are M = ~>.(u.\ I and Mf = Mj, where 
Mj (= 1O~ 3 M ) is a Jupiter mass. Growth timescale tf = 10 5 yr. After t > t{, we set 
M = M f = const. 
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Fig. 3. — Orbital evolution of a swarm of a planetesimals on the b-(e/h) plane. We adopt 
p = 2, Tdamp = 10 4 yrs and U = 10 5 yrs. The planet is fixed at e/h = b = (a p =5 AU). The 
horizontal axis b expresses (a — a p )/h where a is the semimajor axis of planetesimals. Solid 
and dotted lines represent the boundaries of the feeding zone (i.e. Jacobi energy Ej = 0) 
and those at t — 0, respectively. The time evolution of the latter is caused by increase in 
h. The selected number of planetesimals is 1000 in 3.3 AU < a < 8.4 AU at t — 0. The 
numbers of planetesimals are 998(10 3 yr), 992(3 x 10 4 yr), 904(10 5 yr), and 878(3 x 10 5 yr). 
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Fig. 4. — Evolution of the planetesimal accretion rate onto the growing planet as a function 
of the protoplanet mass (M). (a) The results in the gas-free case, the cases of (b) Td amp = 10 6 
yrs, (c) 10 5 yrs and (d) 10 4 yrs. The four lines in the each panel represent the results with 
various gas accretion models (p = 2, 1,0, —2). Initial mass of the protoplanet Mo is set as 
5.67M e . The systems initially consist of 20,000 planetesimals, so the individual planetesimal 
masses correspond to ~ 7.7 x 10 _4 M e . 
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Fig. 5. — Evolution of the planetesimal accretion rate as a function of £ = f///f S cat for 
p — 2, 1, and —2. (a) The results in the gas- free case, the cases of (b) Td amp = 10 6 yrs, (c) 
10 5 yrs and (d) 10 4 yrs. The fitting formula, eq. ff24"]) . is expressed by thick solid lines in the 
plots. 
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Fig. 6. — Evolution of the scaled planetesimal accretion rate as a function of rj = f^/t^amp 
for p = 2, and —2 in the case of Tdamp = 10 4 yrs (left panel) and for Td am p = 10 6 , 10 5 and 
10 4 yrs in the case of p = — 2 (right panel). The fitting formula, eq. (125]) . is expressed by 
thick solid lines. 
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Fig. 7. — Comparison between the numerical simulations and the semi-analytical results. 
The left and right panels show the results for gas accretion with p = 2 and p = 0, respectively. 
The thin and thick curves represent the numerical and semi-analytical results. 
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Fig. 8. — The evolution of cumulative mass of accreted planetesimals as a function of M 
in the case with M = 10M e . The left and right panels show the results for a = 5.2 AU 
and a = 9.55 AU, which correspond to Jupiter and Saturn. Here, R = 2R 1 and fa = 2 are 
assumed, where R\ is the physical radius for mass M and p = lgcm -3 . For other R and fa, 
the accreted mass is multiplied by (R/2Ri) 1 ^ 2 (fd/2). 
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Fig. 9. — The parameter range in which phase 2 can occur, which is expressed by the shaded 
regions. The region above the solid line represents rj > 1. Phase 2 can occur in the regions 
above the dashed line for rj < 1 and in the regions below the dot-dashed line for rj > 1. 



